#
# ==============
# Figure 2: Improper payment rate trends and model fixed effect estimates by fiscal year.
# ==============
#

  AG = data.table(read_dta(file.path(DATALOC,"AgencyLevel.dta")))
  PR = data.table(read_dta(file.path(DATALOC,"ProgramLevel.dta")))

  # Annual unweighted and dollar-weighted mean IP rates.
  FY = AG[FY > 2004,.(Mean_IP=mean(iprate,na.rm=T),ip_dollars=sum(improperpaymentamountm,na.rm=T),outlays=sum(outlaysamountm,na.rm=T)), keyby="FY"]
  FY[, Weighted_mean_IP := 100*ip_dollars/outlays]
  FY[,c("ip_dollars","outlays") := NULL]

  # Estimate TWFE agency level and extract FEs.
  lm_ag = AG[FY >= 2004,lm(iprate ~ agen_hear_all + last_agen_hear + factor(agency) + factor(FY))]
  print(summary(lm_ag))
  fred = confint(lm_ag,level=.95)
  larry = coef(lm_ag)
  FY[,lb95_agen := fred[grepl("FY",rownames(fred)),1]]
  FY[,coef_agen := larry[grepl("FY",names(larry))]]
  FY[,ub95_agen := fred[grepl("FY",rownames(fred)),2]]

  # Estimate TWFE program level and extract FEs.
  lm_pr = PR[FY >= 2004,lm(iprate ~ agen_hear_all + last_agen_hear + factor(programname) + factor(FY))]
  #print(summary(lm_pr))
  moe = confint(lm_pr,level=.95)
  curly = coef(lm_pr)
  FY[,lb95_prog := moe[grepl("FY",rownames(moe)),1]]
  FY[,coef_prog := curly[grepl("FY",names(curly))]]
  FY[,ub95_prog := moe[grepl("FY",rownames(moe)),2]]

  print(FY)

  # Vertical lines and arrows for statutes and EO.
  # Federal government fiscal years run from October 1 of the calendar year prior through September 30 of that year.
  events = data.table(event=c("Obama EO 13520","OMB M-10-13 (EO 13520)","IPERA 2010","IPERIA 2012","Threshold reduced to 1.5%","PIIA 2019"),date=c("2009-11-20","2010-03-22","2010-07-22","2013-01-10","2013-10-01","2020-03-05"),FY=c(2010,2010,2010,2013,2014,2020))

  # Plot. Top panel means, bottom panel FEs.
  # Top panel.
  f = file.path(FIGURELOC,"Figure02-PanelA.pdf")
  pdf(f,width=10,height=6)
  par.old = par(mar=c(3.1,3.1,1.1,0.1),cex.lab=1.1,cex.axis=1.1,mfrow=c(1,1))
  FY[,plot(x=range(FY),y=range(cbind(Mean_IP,Weighted_mean_IP)),type="n",ann=F,axes=F)]
  x_at = seq(2005,2021,2)
  axis(1,at=x_at,labels=x_at);axis(2,las=2);grid(lwd=2)
  events[,abline(v=FY,lwd=2)]
  title(ylab="Improper payment rate",line=2)
  # Lines for mean and weighted mean.
  FY[,lines(x=FY,y=Mean_IP,lwd=4,col=1)]
  FY[,lines(x=FY,y=Weighted_mean_IP,lwd=4,lty=2,col=2)]
  legend("topleft",legend=c("Mean IP rate","Dollar-weighted mean IP"),cex=1.2,col=c(1,2),lty=c(1,2),lwd=4,bg="white")

  # Identify each of the events plotted.
    anArrow <- function(x0,y0,x1,y1,length=.2,label=NULL,...) {
      # Make an arrow on a plot with optional text.
      arrows(x0,y0,x1,y1,length=length)
      if (!is.null(label)) {
        text(x0,y0,label=label,...)
      }
    }
  # FY 2010.
  anArrow(x0=2007.5,x1=2009.9,y0=4.2,y1=3.2,label="Exec Order 13520,\nOMB M-10-13,\nIPERA 2010",pos=2)
  # FY 2013.
  anArrow(x0=2011.5,x1=2012.9,y0=6.4,y1=5.3,label="IPERIA 2012",pos=3)
  # FY 2013.
  anArrow(x0=2015.2,x1=2014.1,y0=7.2,y1=6.2,label="OMB\nreduces\nthreshold\nto 1.5%",pos=4)
  # FY 2020.
  anArrow(x0=2018.8,x1=2019.9,y0=2.7,y1=3.5,label="PIIA 2019",pos=1)
  par(par.old)
  dev.off()


  # Bottom panel.
  f = file.path(FIGURELOC,"Figure02-PanelB.pdf")
  pdf(f,width=10,height=6)
  par.old = par(mar=c(3.1,3.1,1.1,0.1),cex.lab=1.1,cex.axis=1.1,mfrow=c(1,1))
  FY[,plot(x=range(FY),y=range(cbind(lb95_prog,ub95_prog)),type="n",ann=F,axes=F)]
  x_at = seq(2005,2021,2)
  axis(1,at=x_at,labels=x_at);axis(2,las=2);grid(lwd=2)
  #events[,abline(v=FY,lwd=2)]
  title(ylab="Improper payment rate year fixed effect",line=2)
  title(xlab="Fiscal year",line=2)
  # Points.
  FY[,points(x=FY-.1,y=coef_agen,col=3,cex=1.5,pch=19)]
  FY[,points(x=FY+.1,y=coef_prog,col=4,cex=1.5,pch=15)]
  # Segments.
  FY[,segments(x0=FY-.1,y0=lb95_agen,y1=ub95_agen,col=3,lwd=2)]
  FY[,segments(x0=FY+.1,y0=lb95_prog,y1=ub95_prog,col=4,lwd=2)]
  legend("topleft",legend=c("Agency model","Program model"),title="Fiscal year fixed effects",cex=1.2,col=c(3,4),pch=c(19,15),bg="white")

  par(par.old)
  dev.off()